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Abstract. We review our recent work leading to steady-state solutions of the 
semiclassical (Maxwell-Bloch) equations of a laser. These are coupled non-linear partial 
differential equations in space and time which have previously been solved either by 
fully time-dependent numerical simulations or by using major approximations which 
neglect non-linear modal interactions and/or the openness of the laser system. We have 
found a time-independent technique for determining these stationary solutions which 
can treat lasers of arbitrary complexity and degree of openness. Our method has been 
shown to agree with time-dependent numerical solutions to high accuracy and has been 
applied to find the electric field patterns (lasing modes) of random lasers, which lack 
a laser cavity and are so strongly damped that the linear system has no detectable 
resonances. Our work provides a link between an important non-linear wave system 
and the field of quantum/wave chaos in linear systems. 
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1. Introduction: semiclassical laser theory and its limitations 

Lasers are among the most important non-linear systems in modern science and 
technology. They take advantage of the quantum phenomenon of stimulated emission 
to create an amplifier for light, which becomes an oscillator emitting in one or 
several very small frequency windows above a series of thresholds. Shortly after the 
demonstration of the first lasers in the early 1960's (following experiments and theory on 
the maser and proposals for the laser), standard theoretical descriptions were developed 
independently by Haken and Sauermann [1] and Lamb and collaborators [2], [3] which 
neglect the quantum fluctuations of the electromagnetic field, hence these are referred 
to as semiclassical theories. Shortly thereafter a full quantum theory of the laser was 
developed by Haken [1], and Scully and coworkers [5], which builds on the semiclassical 
theories. 

The semiclassical theory describes the electric and magnetic fields using Maxwell's 
equations, which couple the fields to matter through the non-linear polarization of 
the gain medium, which is described by quantum equations of motion. This theory 
describes all of the characteristic non-linear phenomena observed in lasers, such as mode 
competition and selection, bistability, frequency and phase locking. The semiclassical 
laser equations in their most common form describe a gain medium of identical two- 
level "atoms" with energy level spacing hua and relaxation rate 7|| being pumped by 
an external energy source, Dq, contained in a cavity which can be described by a linear 
dielectric function, e(x). This leads to a population inversion of the atoms, -D(x, t) which 
in the presence of an electric field creates a non-linear polarization of the atomic medium, 
P(x, t), which itself is coupled non-linearly to the inversion through the electric field, 
E{'x,t). The electric field and the non-linear polarization are related linearly through 
Maxwell's wave equation, although the polarization is implicitly a non-linear function 
of the electric field in a manner we will demonstrate explicitly below. The induced 
polarization also relaxes at a rate 7^ which is typically much greater than the rate 7|| 
at which the inversion relaxes. This is a key assumption in our approach, as well as in 
almost all of the work in the literature on multi-mode lasing. 

The resulting system of non-linear coupled partial differential equations for the 
three fields E(x, t), P(x, t), L)(x, t) are (c = 1): 

= ^V^i?+ - 4^P+ (1) 
e(x) e(x) 

P+ = - (lUa + 7±)^"' + ^E+D (2) 

in 

b = 71 {D, -D)-^ (e+(P+)* - P+{E+y) . (3) 

Here g is the dipole matrix element of the atoms and the units for the pump are 
chosen so that Dq is equal to the time-independent inversion (uniform in space) of 
the atomic system in the absence of an electric field. The electric field and polarization 
are real functions (vector functions in general, but we assume a geometry where they 
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can be treated as scalars). In writing the equations above we have written these fields 
in the usual manner in terms of their positive and negative frequency components, 
E = + E~ , P = + P~, and then made the rotating wave approximation (RWA) 
in which the coupling terms correspond to driving far away from the atomic transition 
frequency are neglected. This approximation is both very good under general conditions 
and not essential; it just simplifies the equations. We have not made another widely used 
approximation, the slowly-varying envelope approximation. In this approximation the 
fields are factorized as an oscillatory term at the atomic transition frequency multiplied 
by a slowly varying envelope function, whose second derivative is then neglected. We 
have shown recently [6] that this approximation is not very accurate for micro cavities, 
and it provides no significant simplification in our approach. 

Note that if in this non-linear system the electric field and polarization are zero 
at t = then they remain zero and the pump just finds its steady-state value 
i5(x, t) = Dq. This is because the semi classical theory neglects spontaneous emission, in 
which the atom spontaneously relaxes from the upper level through the emission of light, 
generating a fiuctuating electric field. If spontaneous emission is not neglected, then the 
pump will amplify this emission leading to a broad enhancement of light emission near 
the atomic transition frequency. At a series of thresholds for multi-mode lasing much 
narrower lasing lines emerge out of this broad peak of amplified stimulated emission. 
This process is simplified in the semiclassical theory. 

In the semiclassical theory if there is an initial pulse of electric field to start the 
non-linear interactions, then three situations can occur: 

(i) The pulse will decay as t ^ oo and there is no field in the cavity in the presence of 
the pump. In this case the pump rate is below the lasing threshold. 

(ii) There is a transient regime of a fiuctuating field after which the system self-organizes 
to oscillate and emit light at a discrete set of lasing frequencies. When there is only 
one frequency we are in the single-mode regime; otherwise we have multi-mode 
lasing. In the single-mode regime the inversion -D(x, t) will no longer be uniform in 
space, but will be independent of time. In the multi-mode regime the inversion will 
in principle depend on both time and space. The polarization will oscillate with 
the same frequencies as the field. 

(iii) The fields never settle down to a multiperiodic time dependence and quasiperiodic 
or chaotic behaviour occurs. We neglect this possibility in the approach given 
below. Simulations of ([I])-® indicate this does not occur in the parameter range 
we are investigating [6]. It is well-known, however, that chaotic time-dependence 
can occur in multi-mode lasers in general [7]. 

Above the first lasing threshold in the semiclassical theory the inversion becomes 
non-uniform in space and is reduced from its value Dq in the absence of lasing emission. 
This changes the effective gain in the cavity in a space-dependent manner, making it 
more difficult for higher lasing modes to reach threshold. Also, the single mode that is 
lasing saturates as the pump is increased. Thus the modes of a laser have an effectively 
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infinite order non-linear interaction through the gain medium; this effect is called spatial 
hole-burning. Unlike previous theories, our approach is able to treat this effect to all 
orders, in an approximation which is found to be excellent. 

The laser is an open system, and the electric field in a laser is being generated by the 
stimulated emission of light from the atomic medium in the absence of any input light at 
the lasing frequencies. As a consequence the stationary electric field in the laser cavity 
does not conserve energy; the flow of light energy (magnitude of the Poynting vector) is 
increasing along the cavity in the direction of the output (s). Therefore the electric field 
in the cavity has a non-hermitian character and we can show that it is determined by 
a non-hermitian boundary condition; outside the cavity there is no gain medium and 
energy flux is conserved. This non-hermitian character of the actual lasing modes has 
been widely neglected in semiclassical laser theory which has until very recently only 
dealt with hermitian closed cavities in the multi-mode regime [7]. Heuristic arguments 
were used to discuss the output power of the laser and its emission pattern. Modern 
nanofabrication capabilities along with potential applications have led to a wide variety 
of new microlasers with complex cavity structures; examples are photonic bandgap 
defect mode lasers [H], deformed dielectric cavity lasers [9] and random nano-composite 
lasers [lO]. For such lasers the emission pattern and output power is not easily guessed 
from internal mode properties of the closed resonator. Thus we sought a more rigorous 
and predictive framework for analyzing arbitrarily complex and open multi-mode lasers. 
Our approach [6], [11], which we call Ab Initio Self-Consistent (AISC) laser theory is able 
to treat lasers with any degree of openness, and as mentioned earlier, with arbitrarily 
strong multi-mode non-linear interactions. Below we will review this new approach and 
describe a dramatic application to random lasers. 

Random lasers fiUi [12] are aggregates of nanoparticles which have gain, or are 
embedded in a gain medium. These lasers have no mirrors or reflecting boundaries 
to confine light; light generated in the gain medium by stimulated emission simply 
multiply-scatters as it diffuses out to the boundary and escapes from the aggregate. We 
refer to these systems as diffusive random lasers (DRLs). An original motivation for 
studying random lasers was the interest to possibly observe Anderson-localized states of 
the linear wave equation; however such localized-modes random lasers have only been 
demonstrated in one dimension [13]. DRLs correspond to the limit of extreme openness 
where the system in the absence of gain has no long-lived resonances at all and hence 
no detectable natural frequencies to determine the oscillation frequencies of the laser. 
Nonetheless random lasers behave in most respects like conventional multi-mode lasers, 
with sharp thresholds and line-narrowing above threshold. Until our work there had 
been no solution of the semiclassical lasing equations for DRLs above threshold, except 
for numerical simulations of the equations in time and space which provides no physical 
picture of the lasing modes in such systems. Below we show that our AISC laser theory 
is able to find and interpret the lasing modes in such a system, which have been beyond 
conventional treatments of the semiclassical lasing equations. 
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2. Self-consistent steady-state lasing equations 

The starting point of our new formulation is to assume that there exists a steady state 
multiperiodic solution of equations ([I])-© above, for long times after the pump is turned 
on, i.e., we try a solution of the form: 

N N 

i?+(x, t) = ^ ^^(x)e-^*=-*, P+(x, t) = Y^ P^(x)e-^'=-* . (4) 
fj-=i M=i 

Having taken c = 1 we do not distinguish between frequency and wavevector. The 

functions ^'/^(x) are the unknown lasing modes and the real numbers /c^ are the unknown 

lasing frequencies. Unlike standard approaches, here these functions and frequencies are 

completely general and bear no a priori relationship to the normal modes or linear quasi- 

modes of the passive cavity (cavity without gain). 

Note from this ansatz that the non-linear term in ([3]) for the inversion ~ E~^{P~^)* is 

time independent in steady-state for a single- mode solution (A^ = 1), allowing an exactly 

time independent inversion (which does vary in space). However when there are two or 

more modes, the non-linear term in ([3]) will now oscillate at the difference frequencies 

of the modes, ruling out a strictly stationary inversion in the multi-mode case. This 

beating represents the interference of the different modes since they oscillate at different 

frequencies. The other time scale in the dynamics of the inversion is 7|~^; if this time 

scale is long compared to the beat periods of the modes, i.e., if 7|| << AA;^,^ = \kfj_ — k^\, 

\ffi, u, then these time-varying interference terms average to zero. In addition it must 

be assumed that the polarization follows the inversion closely, so 7_l >> 7||. In these 

two limits the inversion may be approximated as stationary to high accuracy. This is 

a standard approximation, often used in multi-mode laser theory and often realized in 

lasers of interest; we will refer to it as the stationary inversion approximation (SIA). 

Once this approximation is made, then ([2]) can be solved to give the polarization in terms 

of the lasing modes ^l/^(x) and the unknown (stationary) inversion, Ds{x.). Substitution 

of this result back into the equation for the inversion gives -Ds(x) in terms of {\E'^(x)} so 

that the polarization in ([1]), which is the source term for the wave equation, is completely 

given in terms of {\E'^(x)}. The time derivatives needed in ([T]) are determined by the 

multiperiodic ansatz (jlj) and one obtains an inhomogeneous differential equation in space 

only for E~^. This nominally inhomogeneous linear differential equation can be formally 

inverted with a Green function (determined by the outgoing wave boundary conditions) 

to yield a self-consistent equation determining all non-zero lasing modes. The details of 

this have been given in P, [TT] ; the resulting equation is: 

v|> = [ , . Do{^)G{^,^-k^)m,{^) 

^ lL-i{k,-k,)kl]v e(x')(l + E.r,|vl/,(x')P)- 

Here ka is the frequency of the atomic transitions, G(x, x'; k^) is the Green function of 
the open cavity, = T^k,^) is the gain profile evaluated at k^, -Do(x) = -Do(l + '^o(x)) 
is the pump, which we have now allowed to vary in space, e(x) = ?7,^(x) is the dielectric 
function of the " cavity" , which can be arbitrarily complex or leaky. This linear dielectric 
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function of the cavity can also be frequency-dependent and complex, but we do not 
consider these cases here. Electric field and pump strength are measured in the natural 
dimensionless units of the problem, Cc = h^'j±'y\\/2g, Dqc = ^7±/4vr/c^5f^. The integral 
in ([5]) extends over the gain region, denoted by V. 

These are the fundamental equations of the AISC laser theory, a set of self-consistent 
non-linear integral equations for the unknown lasing modes and frequencies. We have 
developed an efficient method for solving ([5]) iteratively, which will be described below. 
As expected, these equations exhibit a series of thresholds D^j^ in the pump strength, Dq. 
\i/^j(x) = is always a solution, but above some value of Dq a first non-trivial solution 
appears. It then interacts with itself through the non-linear hole-burning denominator 
in ([5]) and at the same time partially suppresses other solutions which would turn on 
in the absence of the first mode. As the pump increases additional solutions do appear 
in general, with those solutions favored which overlap less in space with the existing 
lasing mode(s). Physically this is because gain is less depleted where the electric field of 
the existing lasing modes is small. Note the presence of the unknown lasing frequencies 
both explicitly and implicitly in ([5]). They are determined by the requirement that the 
non-linear integral operator on the right hand side of ([5]) have real eigenvalues for each 
solution, and the lasing frequency must be allowed to fiow with each iteration of ([5]) to 
maintain this "gauge" condition. More details of this procedure will be given below. 
Despite the fact that this is a non-linear system, there is a natural basis set in which 
to expand the functions \E'^(x), which is determined by the Green function in ([5]). This 
basis set has the property that for a high-Q cavity each lasing mode become equal to a 
single basis function, as will be discussed in the next section. In general ([5]) then can 
be formulated as a non-linear map for the vector of coefficients in this basis set and the 
solutions are the fixed points of this map. 



3. Open cavity boundary conditions and CF states 

The Green function which appears in ([5]), G, is that of the cavity wave equation 

[e(x)-iV2 + G(x, x'; k) = 53(x - x'), (6) 

with the real parameter k = being the unknown lasing frequency of mode /i in 
([5]). The crucial point in the above equation is the nature of the boundary conditions 
on G; at infinity only outgoing waves of frequency /c^ are allowed: VrG'(x, x'; A;) = 
Vr/G(x, x'; /c) = ikG{:K,x.'; k), where Vr is the radial derivative. These are non- 
hermitian boundary conditions implying that the spectral representation of G(x, x'; k) 
is of the form: 

G(x,x';.)^i:^^|^^f^. (T) 

We refer to the functions v'mlx, k) in ([7]) as the constant-fiux (CF) states, as they 
conserve the photon fiux outside the cavity. They satisfy 

- e(x)- Wr„(x, k) = kl{k)ipm{^, k) (8) 
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with the corresponding non-hermitian boundary condition of purely outgoing spherical 
waves at infinity. The non-hermitian boundary conditions imply that the CF 
wavevectors km are always complex, and in fact it can be shown that they always have 
negative imaginary part, corresponding to amplification within the cavity [11]. For all 
the cases considered here it is sufficient to impose the outgoing wave boundary condition 
at the boundary of the cavity. In fact the Green function only correctly calculates the 
electric field inside the cavity. By connecting the solutions to outgoing solutions of the 
free wave equations these solutions determine the electric field everywhere. 

Below we will illustrate the theory first with an application to a ID cavity consisting 
of a perfect mirror at the origin, and a uniform dielectric region of real index n = uq and 
length a, terminated with vacuum out to infinity in the positive x-direction. For this 
special case (see the inset of figure H] left panel) the outgoing wave boundary condition 
is just dx^m{,x)\a = +ikiprn{ci). Notc this differs subtly but importantly from the quasi- 
bound state boundary condition, for which the complex eigenvalue km replaces the real 
wavevector k [H]. Thus the quasi-bound (QB) states, which are often thought of as 
becoming lasing modes when gain is added, are not an appropriate basis, since they have 
complex wavevectors outside the cavity, do not conserve energy flux there, and diverge at 
infinity. This is not true for the CF states, as already noted, which conserve flux outside 
the cavity, but play the role of the linear cavity resonances inside the cavity. Note that 
the CF and QB states never have exactly the same complex eigenvalues, and, unlike the 
QB states, the CF states are a parametric family of basis functions, parameterized by 
the lasing frequency, k. When the cavity is very high-Q, there will be a lasing mode 
near the longest lived cavity resonance (under the gain curve), and in this case there 
will be one CF state with nearly the same complex eigenvalue (see figure [H inset). 

Because the CF states are determined by a non-hermitian boundary condition, 
they are not orthogonal to one another, but instead they are biorthogonal to a set 
of dual functions, (fm{x,k), which satisfy the complex conjugate differential equation 
with the complex conjugate boundary conditions (for the ID case: dxi^mix)\a = 
—ikifmiO')) ■ In general these dual sets of functions satisfy the biorthogonality relations 
[H]: Jjy d:Ki^'^{-x, k)(pni^, k) = 6mn, and are also complete. These relations make it 
possible to expand an arbitrary lasing solution, 

oo 

M/,(x) = Y: aM^) , (9) 

m=l 

SO that each \E'^ is defined by this vector of complex coefficients in the space of 
biorthogonal functions. Because only CF states with frequencies near the center of 
the gain curve contribute to the lasing state, it is possible to truncate the sum in ([9]) 
to a finite number (A^) of components, making \E'^ a finite dimensional vector. In the 
calculations we show here 10 < < 25 is sufficient for convergence (in general higher 
values require larger A^) . Recall that the CF boundary condition depends on the lasing 
frequency, k^, so that v^m(x) = ipmi'^ik^) and in what follows we define ki^ = kmik^). 
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Figure 1. Constant-flux (CF) states for simple edge-emitting ID laser (inset, 
figure 4, left panel), which are the natural basis set in which to express the lasing 
solutions of the AISC equation ([5]). CF states satisfy the linear wave equation 
([8]) with the non-hermitian boundary conditions of only outgoing plane waves 
at the lasing frequency outside the cavity boundary. Hence for this case they 
take the form s'm{nkmx) with complex amplifying wavevector km (3 examples 
shown). The complex CF wavevectors are shown in the inset (red plus signs) 
and compared to the quasibound wavevectors (black open circles) for = 20 
(dashed vertical line). Note that the QB and CF wavevectors nearest to the 
lasing frequency coincide, but otherwise they disagree significantly. 



By substitution of iQ into ([5]) and use of the biorthogonality relations one finds: 

^Dol± (kl/kl) f , (l + 4(xO)^/;*(xO 



(7± - - ka)) {kl - kffn^) Jv e(x')(l + E.r(A;,)|^,(x' 

This is the form of that is employed in our algorithm for finding the lasing modes and 
frequencies. As discussed above, it reduces the problem to finding the complex vector of 
coefficients and the frequency k^ for each lasing mode, which depends non-linearly 
on all the other lasing modes and itself through the infinite order non-linearity evident 
in the denominator of equations ([5]), (fTOj) . 



4. Solution method for AISC equations 

The non-linear amplitude equations in (ITO!) can be written in the compact form 

af" = DoTf'a'' (11) 

where {T^)mn = Tmn{k^; {k^, a'^}) is a non-linear operator acting on the properly 

ij*, a2 , . • • , a^) 

which takes the form: 

,( l + do(xO)(^:,(x^A:)y.„(x^ fc) 
e(x')(l + E.r,|M/,(x')|2) 

where Am{k) = ij±{k'^/kl)/[{'j± — i{k — ka)){k'^ — kl^ik))]. Below some finite value 
of the pump Dq only the trivial solution exists: \E'^ = 0,VyU. As Dq is increased, a 



truncated iV-dimensional vector space of complex amplitudes a'^ = (a^, . . . , a^- 
he form: 

Imn{k,ik,,a }}-Am{k)J^d^ efxOfl + E.rjM/JxOP) ' ^^^^ 
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series of thresholds are reached at which the number of non-trivial solutions increases 
by one. If we denote the threshold for z-mode lasing by D^l^ there exist i solutions 
{/c^,a^} (/i = 1, . . . ,i) for pump parameter Dq such that D^)^ < Dq < dIj^^\ In each 
of these intervals we assume that i solutions to ([5]) exist and find them by a method to 
be described below. 

Very close to the first threshold Dq = D^j^ + e, we may consider the linearized 
operator T^^\k) = T{k; {K, a" = 0}) 

which is obtained by neglecting the Y.u^u\'^u{'^')\'^ term in the denominator of (ITT]) . 
The resulting linear equation associated with f|T3l) has the form, 

r(°)(A;)a^ = (l/Do)a^. (14) 

This equation can in general not be satisfied for a real Dq except at discrete values of 
is a non-hermitian matrix and has complex eigenvalues A„(/c) for general 
values of k. As k is is varied, the eigenvalues Xnik) flow in the complex plane, each one 
crossing the positive real axis at a specific A;„, determined by Im[A„(A; = A;„)] = 
(see figure [2]). The modulus of the eigenvalue defines the "non-interacting" lasing 
threshold corresponding to that eigenvalue, -D^^q = 1/Xn{kn) (these real eigenvalues 
will be denoted by \^^^ = \^^\k = k^))^ the real wavevector k^ is the non- interacting 
lasing frequency, and the eigenvector a" gives the "direction" of the lasing solution in 
the space of CF states. Among these solutions, the smallest d\^q (i.e., the largest of the 
real eigenvalues A^f-*) gives the actual threshold for the first lasing mode; the frequency 
ki is the lasing frequency at threshold and the eigenvector a} defines "direction" of the 
lasing solution at threshold. The "length" of cannot be determined from the linear 
equation but rises continuously from zero at threshold and is determined by the 
non-linear equation ffTTl) infinitesimally above threshold. As noted, the remaining real 
eigenvalues of 'T'^'^\k) define the non-interacting thresholds ior other modes, however the 
actual thresholds of all higher modes will differ substantially from their non-interacting 
values due to the non-linear term in (fTT]) which now comes into play. The actual lasing 
frequencies of higher modes have a relatively weak dependence on Dq and differ little 
from their non-interacting values. Above the threshold D^j^ we solve the non-linear 
equation (fTT!) by an iterative method to be described below. 

Assuming we have the non-linear solution for the first lasing mode, \E^i(x), available 
at each value of Dq (we will denote these discrete values by -Do(O) can construct the 
first interacting threshold matrix as we increment Dq > -D^^'': 

rw(k) - A (k) f ,7,,/ (l + ^o(xO)<^;.(x^fc)y„(x^fc) 

the second largest eigenvalue of which will determine the second interacting threshold, 

(2) 

Dli^ . This procedure can be generalized as additional modes turn on to define the 
second, third, etc. interacting threshold matrices. This procedure gives us a way to 
monitor when (fTTl) has a second, third, etc. non-trivial solution. We diagonalize the 
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Figure 2. Eigenvalue flow for the non-interacting threshold matrix (k)- 
Coloured circles represent the trajectories in the complex plane of 1/X^{k) as 
a function of k for a scan k = {ka — 7_l, ka + 7_l) {kaR = 30, 7^ i? = 1). Shown 
are only a few eigenvalues out of = 16 CF states in the basis. The direction 
of flow as k is increased is indicated by the arrow. Star-shapes mark the initial 
values of 1/X^{k) at k = ka — 'y±- The full red circles mark the values at which 
the trajectories intersect the real axis, each at a different value of A; = k^, 
defining the non- interacting thresholds -D|/^q. 



matrix (llSp and plot DQ(i)\^l^ as we increase Dq in small increments solving the non- 
linear equation ffTTl) . The resulting plot is shown in figure [31 When Do{i)X^^^ reaches 
unity, a new mode has reached threshold. One can confirm that the n*'* threshold matrix 
has n eigenvalues equal to unity until threshold, when the (n + 1)*^ appears. Such a 
diagram is very interesting because it shows the strong effects of mode competition. 
The eigenvalues of the non-interacting threshold matrix are fixed numbers, independent 
of pump, Dq, so that a plot DqX^^^ vs. Dq just gives straight lines intersecting unity 
at the non-interacting thresholds. The interacting threshold matrix depends strongly 
on Dq, so that the interacting eigenvalues will be sub-linear, leading to much higher 
thresholds, and some will even be decreasing with increasing Dq, indicating modes 
which are completely suppressed by mode competition and will never turn on. Thus the 
interacting threshold matrices give us access to the mode competition below threshold, 
as well as the lasing frequencies below threshold, even though the modal amplitudes are 
strictly zero below threshold in the semiclassical theory. 

Using the information obtained by these successive linearizations of the lasing map 
f fTT]) we can begin our iteration process of the non-linear problem with a very accurate 
approximation to the lasing frequencies and the vectors a^' defining each lasing mode. 
This allows the iterative solution to converge in a reasonable time. However there 
is one additional subtlety concerning (fTT]) : this equation is invariant under a global 
phase change a'^ — e^'^a^ . Each member of this continuous family of solutions should 
correspond to a real non-linear eigenvalue 1/Dq for each value of the pump. However 
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Figure 3. Evolution of the thresholds and lasing frequencies as a function of 
Do for the random laser discussed below (see figure O same colour coding is 
used). Left: Evolution of the quantities A^*^Z)o(^) as the pump is increased. 

(i) 

X/j, denotes the real calculated at step i of the discretization of the full 
£)o-range. The lasing threshold corresponds to the line X^Dq = 1, so values 
below that are below threshold at that pump. Once a mode starts to lase its 
corresponding value A^Dq is clamped at A^Z)o/-Doc = 1- The full coloured 
lines represent the modes which start lasing within the calculated range of 
-Do- The dashed lines represent the linear variation X^^^ Do^i) / Dqc of the non- 
interacting thresholds; gray lines indicate modes which do not reach threshold 
in this pump range. Right: Analogous evolution of lasing frequencies with 
pump. The non-lasing modes are drawn in dashed lines, turning into full lines 
as the modes begin to lase. Gray lines represent modes which never lase in the 
calculated range of Dq. Note that frequencies of two lasing modes cannot cross 
(see discussion below), but non-lasing mode frequencies can cross. 



if we iterate our trial solution the phase of the iterate continuously evolves and never 
converges, although the amplitudes of the components do converge. To find a unique 
solution the global phase is fixed (we refer to this as the "gauge condition" ) and the lasing 
frequency is allowed to flow on each iteration to maintain this global phase. This is the 
analog for the non-linear problem of allowing the eigenvalues of the threshold matrices 
to flow to the real axis by tuning frequency eigenvalues discussed above, and depicted in 
figure [3l It ensures that the non-linear eigenvalue is unity, and not a complex number 
of modulus unity. In this manner the interacting lasing frequencies can be found above 
threshold and will differ from the non-interacting or threshold values. In practice, we 
choose the gauge in which we set Im a^j^ = 0, where is the largest CF component 
of the eigenvector a^^ of the non-interacting threshold matrix T^^\k). 

To summarize our solution procedure for (fTTl) : the givens are the dielectric function 
of the resonator, e(x), the atomic frequency, ka and gain width 7_|_. 

• The CF states and wavevectors /c^ are found for a range of k values near ka by 
solving the linear non-hermitian differential equation ([H]). 
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• The first threshold matrix T^^'{k) is constructed from the CF states (and their 
adjoint partners) and frequencies following f|T3|l . k is varied and its largest real 
eigenvalue is found at k = kmax] this defines the first threshold and the lasing 
frequency ki = k^ax] the corresponding eigenvector determines the trial solution 
just above threshold. 

• The pump Dq is increased in small steps and the solution from the previous step 
is used as the trial solution for the next step to solve the non-linear equation (ITTi) 
efficiently by iteration. The lasing frequencies are allowed to flow to maintain the 
gauge condition as discussed above. 

• The solution for the non-linear equation with n lasing modes is used to construct 
the n^^ interacting threshold matrix, T^"\k), which is continuously monitored for 
the appearance of an additional eigenvalue DgX^k) = 1, signaling the threshold for 
mode [n + 1), and the need to increase the number of solutions to the non-linear 
equation f|TT]) . 

Codes based on this algorithmic structure have been developed successfully for 
arbitrary ID cavities, for 2D uniform dielectric cavities of general shape, and for 2D 
disordered cavities embedded in a disk-shaped gain medium. 

5. Solution of ID cavity and comparison to time-dependent solutions 

The advantage of our new approach to semiclassical laser theory is that it gives a time- 
independent theory of the stationary lasing states, which directly calculates the physical 
quantities of interest (thresholds, lasing frequencies, electric flelds inside the cavity and 
outside). The method is computationally much more efficient than fully time-dependent 
simulations, and it provides a physical picture of the lasing solutions as well as analytic 
results near thresholds. Previous time-independent approaches were not quantitative 
and didn't treat the openness of the system fully. Our method is based on one key 
approximation, that of stationary inversion, which holds when 7||/7_l, 7||/AA;^;^ <^ 1, 
where 7|| is the relaxation rate of the inversion and Ak^i, is the frequency spacing of any 
two lasing modes. We tested the accuracy of our approach recently [6] by comparing our 
AISC results for a simple ID edge-emitting laser with exact time-dependent simulations 
of the Maxwell-Bloch (MB) equations for the same system (no assumption of stationary 
inversion). In the latter case the system was simulated until steady-state was reached, 
and then the resulting oscillatory flelds were Fourier-analyzed to yield modal intensities 
and frequencies to compare with the quantities \E'^(x), k^ found with the AISC laser 
theory. 

For the ID edge-emitting laser of length a, uniform index, no, with perfect mirror 
at the origin (see flgure |H left panel) the solutions for the CF states are simply. 




(16) 
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Figure 4. Left: Results of AISC calculation for ID edge-emitting laser with 
perfect mirror at the origin and a dielectric step (n= 1.5— >n = 1) at x = a 
(inset). At this pump level, Dq = 10, three modes are lasing with different 
frequencies and different spatial structures. It is apparent from their growth to 
the loss boundary that these modes do not conserve photon flux and are not the 
solution of a hermitian problem inside the cavity, although they do conserve 
flux outside {x/a > 1). Right: Data points are MB time-dependent solutions 
Fourier-analyzed at stationary output power, solid lines are results of AISC 
calculation requiring 10~^ times the computational effort. Color coding is 
the same as for the left panel. Inset shows lasing frequencies compared using 
the two approaches (diamond vs. solid coloured lines); black curve is the gain 
profile. 



where rim{k) is a /c-dependent normalization constant chosen to give Smn in the 
biorthogonality relation. The complex wavevectors km{k) are found through the 
characteristic equation, 

taninkma) = —i—^ , (17) 
k 

and always have a non-zero negative imaginary part, corresponding to amplification from 
left to right. The lasing solutions are found by the procedure outlined in the section 
above. A single CF state is not an adequate solution for uq = 1.5, but typically only 
contributions from three CF states are needed and the first lasing mode has primary 
contributions from the CF states closest to the atomic frequency, ka- In figure H] we 
show typical results of the AISC laser theory, as well as a comparison of AISC and 
MB solutions. Excellent agreement is found with no free parameters, well into the 
regime of multi-mode lasing in which the inversion is not strictly stationary. This is 
the first demonstration, to our knowledge, of successful quantitative agreement between 
solutions obtained by time-independent methods and the full time-dependent solution 
of MB equations. 

In |6] we varied the parameters 7|[, no and good agreement was found over a wide 
range; in addition, deviations when 7|| — ^ 7_|_ were shown to be amenable to perturbative 
treatment. The time-dependent solutions were full- wave Maxwell-Bloch, because we 



Ab initio self- consistent laser theory and random lasers 



14 



found that the slowly- varying envelope approximation was not very accurate [B]; hence 
these simulations were required to cover the large difference in time scales between rapid 
oscillations at and the slow equilibration at 7|~^, making them very time-consuming. 
In contrast, our theory is more accurate as 7|^^ increases compared to other dynamical 
time scales and requires no more computational effort. Therefore we believe our method 
may be the only practical numerical technique for solving the MB equations for more 
complex and more realistic three-dimensional laser structures. In addition, our method 
is infinite order in the non-linearity and we found [6] that a 3rd order treatment of the 
non-linear interactions failed badly. 




X / a X / a 



Figure 5. Left: The spatial profile of a lasing mode just above threshold (blue 
curve) and well above threshold (red curve). Non-linear self-interaction has 
mixed in higher CF states, increasing the amplification rate for this mode. 
The gain width 7^ = 0.5 was chosen to maintain single-mode lasing over wide 
pump range. Right: Effect of spatial hole burning interaction on higher modes: 
Red curve is the stationary inversion well above threshold in the multi-mode 
regime for the ID laser (inset). Blue curve is the second lasing mode above 
threshold, compared to its profile (green) in the absence of interactions. Note 
the reduction in spatial oscillations due to interactions. The results are scaled 
so that the maximum value is unity. 

As the lasing solutions evolve above threshold, they interact with themselves and 
with other modes in a space-dependent manner. The electric fields burn a hole in the 
inversion which acts as a (complex) index variation for the other modes, and for itself. As 
a consequence the "shape" of the modes continues to change substantially with pump. 
This effect is neglected in conventional approaches which assume each lasing mode is 
a single fixed cavity mode, with a varying overall scale factor, figure [5] illustrates this 
effect. 

6. Diffusive random lasers 

The diffusive random laser (DRL) is perhaps the most challenging new system of 
interest for fundamental laser physics. In its most basic realization, it consists of a 
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random aggregate of particles which scatter hght and have gain or are embedded in 
a background medium with gain [12, 15-23, 25-27]. The system has no traditional 
resonator to trap the light; photons generated by stimulated emission are only very 
minimally "confined" by multiple scattering as it diffuses out of the medium. The 
diffusive escape is so rapid that the medium exhibits no isolated resonances at all in the 
absence of gain. Standard approaches to multi-mode laser theory assume each lasing 
mode is described spatially by a linear cavity mode, usually a closed cavity mode, and 
thus were unsuited to deal with such an extremely leaky cavity. Despite the lack of sharp 
resonances, the laser emission from DRLs was observed to have the essential properties 
of conventional lasers: the appearance of coherent emission with line-narrowing above 
a series of thresholds, and uncorrelated photon statistics above threshold indicative of 
gain saturation [121 US [HI El]. Earlier numerical work (for a recent review see [TO] ) 
has suggested: 1) That the lasing modes are coherent, and not the result of incoherent 
amplification as originally supposed. 2) That these modes exist outside of the strong 
disorder regime, where high-Q linear modes arise due to Anderson localization. These 
observations raise the obvious questions: what are the nature of the lasing modes in 
DRLs? What determines the lasing frequencies, which are not determined by long-lived 
cavity resonances? What is the effect of modal interactions in such complex and low 
finesse cavities? 




Figure 6. Three dimensional rendering of the electric field from AISC numerical 
calculations in the random laser discussed below (see figure [9|) , in which the 
yellow spheres represent the nano-particles in a cylindrically symmetric gain 
medium. The electric field, where color and height indicate intensity, represents 
the steady state solution of the Maxwell-Bloch equations at Dq/Dqc = 123.5. 
The system is illuminated uniformly by incoherent light, shown in this figure 
as coming from above. 
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The AISC theory is the first time-independent approach which can solve the 
semiclassical lasing equations for a multi-mode DRL; a typical example for a calculation 
of the multi-mode electric field both inside the "cavity" and outside is rendered in three- 
dimensions in figure [61 Our first results [28] provide a framework for thinking about the 
lasing modes in DRLs. In [28] we modeled a 2D DRL by a disk-shaped gain medium 
which contains a random distribution of weakly scattering sub- wavelength particles. 
Such a system has strongly-overlapping linear resonances; in the language of resonator 
theory, its finesse, /, (the ratio of the resonance spacing on the real axis to the typical 
distance of a resonance from the real axis) is much less than unity. The AISC laser theory 
indicates that in such a situation the lasing modes will be a superposition of many CF 
states with roughly equal weight; an example of this is shown in figure [71 Thus, as 
compared to the edge-emitting laser, a larger number of modes N ^ 1/f have to be 
kept in the expansion ([9]). The CF states of a DRL around a given ka have comparable 
decay rates = Im [k^] and a similar statistical distribution as the resonances [28] . 




CF3: 26.9% CF6: 17.8"o CF 12:8.4% CF 11: 7.2% 



Figure 7. Top left: false colour intensity plot of typical mode of DRL; example 
is mode 3 at Dq/Dqc = 123 in figure [H Bottom: four CF states which make the 
largest contribution to the lasing mode. White line is boundary of disk-shaped 
gain region. Note that while the CF states fluctuate in space, their variance 
is largest at the cavity boundary. This is due to their non-hermitian nature 
and is discussed below. Top right: schematic illustrating composite nature 
of DRL lasing modes. One mode with a real frequency has contributions 
from roughly 1// CF states with complex wavevector of the type shown in the 
bottom panels. 

Large DRLs with kaR 3> 1 (-R is the linear scale of the random aggregate) are 
typically highly multi-mode laser oscillators at moderate pumping strengths, due to the 
existence of many CF states with similar "decay rates". Km- As a consequence, modal 
interactions through spatial hole burning are extremely strong and have to be taken 
into account correctly. This manifests itself first and foremost in the large shifts of 
the interacting thresholds from their non-interacting values calculated from 'T'^'^\k), as 
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shown in figures [S] and [3] (above) . 
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Figure 8. The shift of the exact thresholds from their non- interacting values. 
The bottom end of the blue lines denote the thresholds determined from the 
eigenvalue flow of T^^^ (/c) while the top of the lines mark the exact thresholds 
calculated by solving ([5]). The red dot at mode number 1 indicates that this 
threshold does not shift. 
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The frequency of a lasing line for a DRL is collectively determined. In [2H] we 
derived a relation of the form: 

K ^ + kf (18) 

where k^^^ satisfies the standard relation for the frequency of a single-mode laser [7j, 
which is a weighted average of the cavity and atomic transition frequencies (we take 
the largest CF component at threshold as the cavity frequency), k^^^ is a collective 
contribution due to all the other CF components which has no analog in conventional 
lasers and involves a sum over the off-diagonal elements of T;5in- When k^R ^ 1 the 
collective term dominates and we expect lasing frequencies to be distributed randomly 
(although correlated) under the gain curve. This is roughly what was found (see 
figure [9]) . 

The evolution of intensities in a DRL as a function of pumping strength Dq (figure [H 
right panel) is very different from that of a conventional edge-emitting laser (figure |4]). 
The intensity variation displays a complex non-monotonic behaviour due to strong 
modal interactions. This is simultaneously refiected in the evolution of the lasing 
frequencies, illustrated in the inset (right panel) of figure |9] above. When two lasing 
frequencies approach one another, instead of repelling, one of the frequencies disappears 
(e.g., the "black" frequency in figure [9]), because the corresponding mode is driven to 
zero, an effect which can't happen in a linear system. Recall that these frequencies are 
not the eigenvalues of a random matrix but rather the parameter values k^ which make 
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Figure 9. Left: Frequencies of multi-mode DRL well above threshold (coloured 
lines). Filled circles are complex frequencies of the CF states contributing 
to lasing modes; coloured circles denote state with largest contribution to 
each mode, but in general each mode is a superposition of many CF states. 
Solid black curve is the gain profile. Inset: DRL model with a random 
configuration of scatterers in uniform disk of gain medium. Right: Modal 
intensities vs. pump for same DRL. Note complex non-monotonic behaviour 
due to modal interactions. Inset: lasing frequencies vs. pump above (solid) and 
below (dashed) threshold. 



the random operator Tf^^ have a real eigenvalue. It is possible to analyze this condition 
in terms of the interacting threshold matrices defined above. In order for frequency 
degeneracy to occur this random complex matrix would have to have two degenerate 
real eigenvalues at the same value of A;^. We argue [28] that level repulsion in the 
complex plane prevents this from happening. However, instead of leading to frequency 
repulsion, this effect leads to frequency locking, i.e., the two interacting modes merge 
into one. We can confirm this effect by noting that after one modes "dies", the CF 
decomposition of the surviving mode has large components associated with the mode 
which was driven to zero. So instead of the "exchange of identity" familiar from linear 
level repulsion, we find "merger of identity". In a time- dependent treatment, this kind 
of cooperative frequency locking is well-known, and we studied it in a chaotic cavity 
laser in [29]. Interestingly, our results show that it is even possible for a mode to "die" 
and then to "reincarnate" itself when its frequency moves far enough away from that of 
the dominant mode (see the mode represented by purple line in figure [9]). 

Recent experiments [30] have shown DRL frequencies to be equally spaced and 
have compared them to the eigenvalue distribution of the gaussian orthogonal ensemble 
of random matrices, which our calculations suggest may not be appropriate, since 
the frequencies are not themselves the eigenvalues of a random matrix. The same 
experiments have shown that for a fixed impurity configuration the frequencies are 
rather stable between different pump pulses, while the modal intensities vary widely. We 
expect this behaviour due to the modal interactions. The frequencies only vary slowly 
as a function of pump, whereas the intensities vary strongly as seen in figure [9] above. 
Thus the exact pump strength or spatial profile of the pump would not be expected 
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to change the frequencies much. However, the interactions can ampUfy the effect of 
any small change in lasing threshold of one mode vs. another. In |28j we found that 
a 30% random variation in the pump profile (for fixed impurity configuration) hardly 
changed the lasing frequencies at all well above threshold, but changed the relative 
modal intensities by factors of order unity (see figure fTOl below) . 




Lasing frequency kR r/R 

Figure 10. Left: Intensity and frequency fluctuations in a DRL. For fixed 
impurity configuration we compare lasing frequencies and intensities in the 
multimode regime for uniform pump (solid lines) and uniform pump with 30% 
zero mean spatial noise added (dashed lines). Right: Field distribution of a 
DRL. Inset: False colour plot of electric field intensity of the DRL of figured] 
(seven modes lasing) at pump Dq/Dqc = 123.5 (white circle is boundary of 
gain medium). Note brightest regions appear at the edge of the gain medium. 
Main plot: Radial intensity of CF states contributing to the lasing modes 
averaged over 400 disorder configurations. There is a large non-random increase 
of intensity with r. 

Finally, we consider the spatial variation of the electric field in DRLs (figures [71 
[TOi) . The false colour representation of the multi-mode electric field in the laser has a 
striking property: it is consistently brighter at the edge of the disk than at its center, 
even though the gain is uniform and there are no special high-Q modes localized near 
the edge. To demonstrate that this effect is not a statistical fluctuation associated 
with this particular disorder configuration we have averaged the behavior of the entire 
basis set of CF states over disorder configurations. The result is an average growth of 
intensity towards the boundary. The origin of this is the complex amplifying character 
of the CF states, which lend the same character to the lasing modes which they form 
by superposition (compare to figures [1] and H]) . The non-hermitian boundary condition 
causes "growth to loss", i.e., a tendency for the electric field to be higher near the 
loss boundary where it must compensate for photons leaking out in steady-state. The 
leakier the laser cavity, the stronger this effect is, so it is particularly pronounced in 
DRLs which have finesse much less than unity. 

Note that this effect means that the electric field fluctuations in DRLs will differ 
substantially from the random matrix/quantum chaos fluctuations of linear cavity 
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modes[2I], first because each mode is a superposition of pseudo-random CF states 
and second because these CF states themselves are not uniform on average. We are 
currently studying both the statistical distribution for lasing frequencies and for the 
spatial characteristics of the lasing modes for the DRL. Another system of interest in this 
regard are the dielectric cavity lasers with chaotic ray dynamics [9]. Here the complexity 
just arises from the shape of the boundary, which need not be highly transmissive, so 
that the statistical properties of the lasing modes should be close to those of single CF 
states. 

7. Summary and conclusions 

A laser is an open non-linear system which self-organizes to exhibit multiperiodic 
behavior. In conventional approaches the lasing modes are described by linear cavity 
resonances and frequencies. We have recently found a method to determine the steady- 
state of a laser which does not start from the linear cavity properties, but allows 
the laser to self-consistently "find" the steady-state lasing modes and frequencies. 
These quantities are determined by a set of infinite order non-linear integral equations, 
for which we have developed an efficient, iterative solution algorithm. This time- 
independent self-consistent method has been tested against exact numerical solution 
of the Maxwell-Bloch (MB) lasing equations for a ID cavity and has been found to 
provide highly accurate solutions over a wide and physically relevant parameter range. 
In our time-independent approach the lasing modes are in general a superposition of 
non-hermitian linear cavity states, termed constant-fiux (CF) states. The more open 
the laser cavity, the more terms contribute to this superposition. An extreme example 
is the diffusive random laser (DRL), which is so leaky that is has no detectable linear 
cavity resonances. We have shown that our method can find the lasing modes of such 
a system, and reveals the effects of strong modal interactions. The theory presented 
here is "ab initio" in the sense that it generates all properties of the lasing states from 
knowledge of the dielectric function of the host medium and two parameters of the gain 
medium (the frequency of the lasing transition and its relaxation rate). We refer to it as 
ab initio self-consistent (AISC) laser theory. The method should be applicable to any 
novel laser-cavity system, and if generalized to treat the gain medium more realistically, 
may allow accurate simulation of lasers for design and optimization purposes. 

The linear CF states provide a link between the problem of quantum/wave chaos 
and the steady-state lasing properties of disordered or wave-chaotic resonators. We 
hope to combine these linear methods along with our non-linear self-consistent equation 
to provide a statistical description of random or wave-chaotic lasers. It will also be 
important to incorporate this full non-linear treatment into a quantum theory of open 
laser cavities [32], which at present remains a challenging and only partially solved 
problem. 
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